% Drude model implementation of 1D and 2D DNG.
\chapter{FDTD Implementation of Metamaterials}
\section{Limitations of Standard FDTD Algorithm}
Standard FDTD\index{FDTD} algorithm cannot cater for negative values of permittivity or permeability. This is because of the Courant stability criterion\index{Courant stability criterion}. As soon as the permeability or permittivity becomes less than unity the algorithm will not remain stable. A metamaterial\index{metamaterial} object can be modelled as a dispersive substance using either the Lorentz\index{dispersive model!Lorentz} or Drude dispersive models\index{dispersive model!Drude}. These models can yield negative values of permittivity (or permeability) for certain frequency ranges~\cite{NumericalFDTD-Sibel}. Using these dispersive models, FDTD update equations\index{FDTD!update equations} are modified and permittivity and permeabilities are replaced with terms dependent on frequency of operation ($\omega$).
\section{Drude Dispersive Model}
In ideal conditions the permittivity (and permeability) of a material remain constant for any frequency and throughout the structure of that material. Speed of electromagnetic waves in such a medium remain constant if frequency changes. Additionally, there is no loss in energy as the waves pass through the medium.

In reality, such a material does not exist. Speed of EM waves varies with frequency of operation. Also, there is a loss associated with the material. A material is dispersive\index{dispersive material} if its permittivity or permeability is dependent on frequency~\cite[Ch. 10]{JBSchneiderUFDTD}.

This section follows the treatment of~\cite[Ch. 10, 289--290]{JBSchneiderUFDTD} for deriving permittivity relationship given by Drude model\index{dispersive model!Drude}. The Drude model is described by the following second order differential equation that relates the net force on charges moving under the influence of an electric field and facing an impeding force due to collisions\index{collision} with material.
\begin{equation}
\centering
M\dfrac{d^2 \textbf{x}}{dt^2} = Q \textbf{E}(t) - Mg\dfrac{d \textbf{x}}{dt}
\label{Drude-Model-Second-Order-DE}
\end{equation}
The left side of this equation gives mass times acceleration or the net force on charge. Here $M$ is the mass of charge, $Q$ the amount of charge, $\textbf{E}(t)$ is electric field that may vary with time and $g$ is the damping coefficient. Equation~\ref{Drude-Model-Second-Order-DE} can be converted to frequency domain using the relationships $\partial ^2/\partial t^2 \rightarrow (j\omega)^2$ and $\partial/\partial t \rightarrow j\omega$, obtaining
\begin{equation}
\centering
M(j\omega)^2\hat{\textbf{x}}(\omega)+Mg(j\omega)\hat{\textbf{x}}(\omega)=Q\hat{\textbf{E}}(\omega).
\label{Drude-Model-Second-Order-DE-Frequency-Domain}
\end{equation}
The polarization vector $\textbf{P}$ is given by
\begin{equation}
\centering
\hat{\textbf{P}}=NQ\hat{\textbf{x}}.
\label{Polarization}
\end{equation}
Here, $N$ is the number of dipoles per unit volume. By eliminating $\hat{\textbf{x}}$ from equations~\ref{Drude-Model-Second-Order-DE-Frequency-Domain} and~\ref{Polarization} and rearranging terms polarization can be expressed as
\begin{equation}
\centering
\hat{\textbf{P}}(\omega)=-\epsilon_0\dfrac{\frac{NQ^2}{M\epsilon_0}}{\omega^2-jg\omega}\hat{\textbf{E}}(\omega).
\label{Polarization-Frequency-Domain}
\end{equation}
By letting $\omega^2_p=NQ^2/M\epsilon_0$, the electric susceptibility is obtained as
\begin{equation}
\centering
\hat{\chi_e}(\omega)=-\dfrac{\omega^2_p}{\omega^2-jg\omega}.
\label{Electric-Susceptibility-Drude}
\end{equation}
The relative permittivity in Drude model\index{Drude model $\epsilon_r(\omega)$} is then
\begin{equation}
\centering
\hat{\epsilon_r}(\omega)=\epsilon_\infty-\dfrac{\omega^2_p}{\omega^2-jg\omega}.
\label{er-Drude}
\end{equation}
Setting $g=0$ and $\epsilon_\infty=1$, relative permittivity comes out to be negative for $\omega/\omega_p > 1$ (figure~\ref{DrudeModel_er}). Thus, Drude model can be effectively used to model metamaterials\index{metamaterial} with permittivity or permeability less than one by incorporating it into FDTD update equations.
\begin{figure}[H]
\centering
\includegraphics[scale=0.8, trim=4cm 8.5cm 4cm 8.5cm, clip]{Figures/FigCh03_DrudeModel_er.pdf}
\caption{$\epsilon_r$ plotted against $\omega/\omega_p$ for $\epsilon_\infty=1$ and $g=0$}
\label{DrudeModel_er}
\end{figure}
\section{FDTD Equations in 1D}
Here, the approach outlined in~\cite{Radial-Zhao} for derivation of FDTD auxiliary update equations is followed. Magnetic field $\textbf{H}$ and magnetic flux density $\textbf{B}$ are related by the constitutive parameter $\mu$, expressed as
\begin{equation}
\centering
\textbf{B} = \mu \textbf{H}.
\label{B-mu-H}
\end{equation}
Where, $\mu = \mu_0 \mu_r$. For a lossless, homogeneous, anisotropic and dispersion--less material, relative permeability or $\mu_r$ is a scalar constant. For dispersive material $\mu_r$ is a function of frequency. Similar to equation~\ref{er-Drude}, relative permeability for Drude model\index{Drude model $\mu_r(\omega)$} is given by
\begin{equation}
\begin{split}
\hat{\mu_r}(\omega)&=\mu_\infty-\dfrac{\omega^2_p}{\omega^2-jg\omega}\\
&=\dfrac{\mu_\infty(\omega^2-jg\omega)-\omega^2_p}{\omega^2-jg\omega}.
\end{split}
\label{ur-Drude}
\end{equation}
From equation~\ref{B-mu-H} and~\ref{ur-Drude},
\begin{equation}
\centering
\nonumber \textbf{B} \left( \dfrac{\omega^2-jg\omega}{\mu_\infty(\omega^2-jg\omega)-\omega^2_p} \right) =\textbf{H}
\label{B-muomega-H}
\end{equation}
and
\begin{equation}
\centering
\omega^2\textbf{B}-g(j\omega)\textbf{B} = \mu_\infty\omega^2\textbf{H}-\omega^2_p\textbf{H}-\mu_\infty g(j\omega)\textbf{H}.
\label{B-muomega-H-frequency-domain}
\end{equation}
Frequency domain quantities can be converted to time--domain using the relationships $j\omega \rightarrow \partial/\partial t$ and $\omega^2 \rightarrow - \partial^2/\partial t^2$. Moreover, fields multiplying with $\omega^2_p$ are averaged in time. Second--order difference scheme is used, both for single and double derivatives to keep all the terms in accordance with the second--order nature of whole expression. That would result in easier implementation. Assuming anisotropic nature of material, only $H_y$ and $B_y$ can be considered. Also, $\omega_{pm}$ and $\omega_{pe}$ are used for magnetic and electric field expressions instead of $\omega_p$, whereas, $\gamma_m$ and $\gamma_e$ are used for magnetic and electric collision frequencies\index{collision!frequency ($\gamma$)}, respectively.
\begin{equation}
\begin{split}
\left(\dfrac{B^{n+1}_y-2B^n_y+B^{n-1}_y}{(\Delta t)^2}\right)&+\gamma_m\left(\dfrac{B^{n+1}_y-B^{n-1}_y}{2\Delta t}\right) =\\
&\mu_0\mu_\infty\left(\dfrac{H^{n+1}_y-2H^n_y+H^{n-1}_y}{(\Delta t)^2}\right)\\
&+\mu_0\omega^2_{pm}\left(\dfrac{H^{n+1}_y+2H^n_y+H^{n-1}_y}{4}\right)\\
&+\mu_0\mu_\infty \gamma_m\left(\dfrac{H^{n+1}_y-H^{n-1}_y}{2\Delta t}\right).
\end{split}
\label{2nd-order-B-H}
\end{equation}
By separating $H^{n+1}_y$, the final form of equation is
\begin{eqnarray}
\nonumber H^{n+1}_y &=& a_m\left(B^{n+1}_y-2B^n_y+B^{n-1}_y\right)+b_m\left(B^{n+1}_y-B^{n-1}_y\right)\\
&&+c_m\left(2H^n_y-H^{n-1}_y\right)+d_m\left(2H^n_y+H^{n-1}_y\right)+e_m H^{n-1}_y.
\label{2nd-order-B-H-final-form}
\end{eqnarray}
Where,
\begin{eqnarray}
\nonumber & a_m=\dfrac{4}{4\mu_0\mu_\infty+\mu_0\omega^2_{pm}\left(\Delta t\right)^2+\mu_0\mu_\infty \gamma_m \left(2\Delta t\right)},\\
\nonumber & b_m=\dfrac{\gamma_m\left(2\Delta t\right)}{4\mu_0\mu_\infty+\mu_0\omega^2_{pm}\left(\Delta t\right)^2+\mu_0\mu_\infty \gamma_m \left(2\Delta t\right)},\\
\nonumber & c_m=\dfrac{4\mu_0\mu_\infty}{4\mu_0\mu_\infty+\mu_0\omega^2_{pm}\left(\Delta t\right)^2+\mu_0\mu_\infty \gamma_m \left(2\Delta t\right)},\\
\nonumber & d_m=\dfrac{-\mu_0\omega^2_{pm}\left(\Delta t\right)^2}{4\mu_0\mu_\infty+\mu_0\omega^2_{pm}\left(\Delta t\right)^2+\mu_0\mu_\infty \gamma_m \left(2\Delta t\right)},\\
\nonumber & e_m=\dfrac{\mu_0\mu_\infty \gamma_m\left(2\Delta t\right)}{4\mu_0\mu_\infty+\mu_0\omega^2_{pm}\left(\Delta t\right)^2+\mu_0\mu_\infty \gamma_m \left(2\Delta t\right)}.
\label{1D-Drude-H-scalars}
\end{eqnarray}
A similar derivation can be carried out to obtain the auxiliary update equation\index{FDTD!auxiliary update equation} for electric field given by
\begin{eqnarray}
\nonumber E^{n+1}_x &=& a_e\left(D^{n+1}_x-2D^n_x+D^{n-1}_x\right)+b_e\left(D^{n+1}_x-D^{n-1}_x\right)\\
&&+c_e\left(2E^n_x-E^{n-1}_x\right)+d_e\left(2E^n_x+E^{n-1}_x\right)+e_e E^{n-1}_x.
\label{2nd-order-D-E-final-form}
\end{eqnarray}
Where
\begin{eqnarray}
\nonumber & a_e=\dfrac{4}{4\epsilon_0\epsilon_\infty+\epsilon_0\omega^2_{pe}\left(\Delta t\right)^2+\epsilon_0\epsilon_\infty \gamma_e \left(2\Delta t\right)},\\
\nonumber & b_e=\dfrac{\gamma_e\left(2\Delta t\right)}{4\epsilon_0\epsilon_\infty+\epsilon_0\omega^2_{pe}\left(\Delta t\right)^2+\epsilon_0\epsilon_\infty \gamma_e \left(2\Delta t\right)},\\
\nonumber & c_e=\dfrac{4\epsilon_0\epsilon_\infty}{4\epsilon_0\epsilon_\infty+\epsilon_0\omega^2_{pe}\left(\Delta t\right)^2+\epsilon_0\epsilon_\infty \gamma_e \left(2\Delta t\right)},\\
\nonumber & d_e=\dfrac{-\epsilon_0\omega^2_{pe}\left(\Delta t\right)^2}{4\epsilon_0\epsilon_\infty+\epsilon_0\omega^2_{pe}\left(\Delta t\right)^2+\epsilon_0\epsilon_\infty \gamma_e \left(2\Delta t\right)},\\
\nonumber & e_e=\dfrac{\epsilon_0\epsilon_\infty \gamma_e\left(2\Delta t\right)}{4\epsilon_0\epsilon_\infty+\epsilon_0\omega^2_{pe}\left(\Delta t\right)^2+\epsilon_0\epsilon_\infty \gamma_e \left(2\Delta t\right)}.
\label{1D-Drude-E-scalars}
\end{eqnarray}
For a wave propagating in $z$ direction, FDTD update equations\index{FDTD!update equations} are
\begin{equation}
B^{n+1}_y(k)=B^n_y(k)+\dfrac{\Delta t}{\Delta z}\left(E^n_x(k)-E^n_x(k+1)\right)
\label{1D-B-Update-Equation}
\end{equation}
and
\begin{equation}
D^{n+1}_x(k)=D^n_x(k)+\dfrac{\Delta t}{\Delta z}\left(H^{n+1}_y(k-1)-H^{n+1}_y(k)\right).
\label{1D-D-Update-Equation}
\end{equation}
Equations \ref{1D-B-Update-Equation} and \ref{1D-D-Update-Equation} drive the FDTD algorithm which give future values of $By$ and $Dx$ from past fields. Equations \ref{2nd-order-B-H-final-form} and \ref{2nd-order-D-E-final-form} are auxiliary equations which give future fields $H_y$ and $E_x$ at $n+1$.
\section{Simulation of 1D DNG Slab}
\subsection{Problem Specification}
An electromagnetic wave travelling in $z$ direction is incident on a slab with negative values of permittivity and permeability (DNG)\index{DNG} at the frequency of operation \cite{DNG-Ehud-Ziol}. Sinusoidal wave\index{source!sinusoidal}, Gaussian pulse\index{source!Gaussian pulse} and Ricker wavelet\index{source!Ricker wavelet} are used as sources. Transmission\index{coefficient!transmission ($\tau$)} and reflection\index{coefficient!reflection ($\Gamma$)} coefficients are calculated at the air--slab interface. Refractive index of slab\index{refractive index ($n$)} for a range of frequencies is also computed. By varying parameters, the simulation can be scaled to any desired frequency or wavelength. Matlab\index{Matlab} code of this simulation is given in appendix \ref{App:1DFDTDDNGMatlab}.
\subsection{Implementation Details}
There are three main phases in the implementation:
\begin{enumerate}
\item Initialisation
\item Simulation
\item Post-processing
\end{enumerate}
Figure \ref{1D-DNG-Algorithm} depicts the flow of program.
\subsection{Initialisation}
In this phase, simulation parameters are specified. Simulation parameters are total number of spatial steps, time steps, wavelength and frequency of operation etc. Field arrays are then allocated and initialised. Any additional fields that are required for post--processing\index{post--processing} are also initialised. Arrays containing $\gamma_m$, $\gamma_e$, $\omega^2_{pm}$ and $\omega^2_{pe}$ must be initialised with values specified for each location in the domain.
\subsection{Simulation}
A dry run\index{dry run} must be carried out without any obstacle if reflection coefficient is to be calculated in order to record the incident field where obstacle is to be placed. The actual simulation consists of five steps:
\begin{enumerate}
\item Compute future value of $B_y$ from past values of $B_y$ and $E_x$ (equation \ref{1D-B-Update-Equation}).
\item Compute future value of $H_y$ from $B_y$ calculated in step 1 (equation \ref{2nd-order-B-H-final-form}).
\item Compute future value of $D_x$ from past value of $D_x$ and $H_y$ in step 2 (equation \ref{1D-D-Update-Equation}).
\item Compute future value of $E_x$ from $D_x$ calculated in step 3 (equation \ref{2nd-order-D-E-final-form}).
\item Update additive source at specified source location for current time step.
\end{enumerate}
Electric field snapshots are saved at specified intervals during simulation and played--back as a movie after simulation ends. 
\subsection{Post--processing}
After the simulation ends, reflection and transmission coefficients are calculated at the air--slab interface. Refractive index is calculated in the slab. A time--domain movie of the simulation is played that shows incident wave passing through the slab and undergoing dispersion.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}
	% Initialisation.
	\draw (5cm,0.1cm) node {\textsf{Initialisation}};
	% Large dashed rectangle.
	\draw[line width=1.2pt, rounded corners, dashed] (0cm,-0.2cm) rectangle (10cm,-3.4cm);
	% Simulation parameters.
	\draw[line width=1pt, rounded corners] (0.5cm,-0.4cm) rectangle (9.5cm,-1.2cm);
	\draw (5cm,-0.8cm) node {\textsf{Initialise Simulation Parameters}};
	% Allocate data arrays.
	\draw[line width=1pt, rounded corners] (0.5cm,-1.4cm) rectangle (9.5cm,-2.2cm);
	\draw (5cm,-1.8cm) node {\textsf{Allocate Data Arrays}};
	% Initialise slab parameter arrays.
	\draw[line width=1pt, rounded corners] (0.5cm,-2.4cm) rectangle (9.5cm,-3.2cm);
	\draw (5cm,-2.8cm) node {\textsf{Initialise Slab Parameter Arrays}};
	\draw[line width=1pt, ->, >=stealth] (5cm,-3.6cm) -- (5cm,-4.5cm);

	% Simulation.
	\draw (5cm,-4.9cm) node {\textsf{Simulation}};
	% Large dashed rectangle.
	\draw[line width=1.2pt, rounded corners, dashed] (0cm,-5.2cm) rectangle (10cm,-12.5cm);
	% Dry run without obstacle.
	\draw[line width=1pt, rounded corners] (0.5cm,-5.4cm) rectangle (9.5cm,-6.2cm);
	\draw (5cm,-5.8cm) node {\textsf{Dry Run Without Obstacle}};
	% Actual Simulation.
	\draw[line width=1pt, rounded corners, dashed] (0.5cm,-6.7cm) rectangle (9.5cm,-12.3cm);
	\draw (5cm,-6.5cm) node {\textsf{Actual Simulation}};
	% Compute H fields.
	\draw[line width=1pt, rounded corners] (1.0cm,-6.9cm) rectangle (8.0cm,-7.7cm);
	\draw (4.5cm,-7.3cm) node {\textsf{Compute $B_y$ and $H_y$}};
	% Compute E fields.
	\draw[line width=1pt, rounded corners] (1.0cm,-7.9cm) rectangle (8.0cm,-8.7cm);
	\draw (4.5cm,-8.3cm) node {\textsf{Compute $D_x$ and $E_x$}};
	% Apply source.
	\draw[line width=1pt, rounded corners] (1.0cm,-8.9cm) rectangle (8.0cm,-9.7cm);
	\draw (4.5cm,-9.3cm) node {\textsf{Update Source}};
	\draw[line width=1pt, ->, >=stealth] (4.5cm,-9.7cm) -- (4.5cm,-10.2cm);
	% Diamond.
	\draw[line width=1pt] (4.5cm,-10.2cm) -- (2.5cm,-11cm) -- (4.5cm,-11.8cm) -- (6.5cm,-10.9cm) -- cycle;
	\draw (4.5cm,-11cm) node {\textsf{q = MaxTime?}};
	% No?
	\draw[line width=1pt, ->, >=stealth] (6.5cm,-10.9cm) -- (8.8cm,-10.9cm) -- (8.8cm,-7.2cm) -- (8.0cm,-7.2cm);
	\draw (6.8cm,-10.7cm) node {\textsf{No}};
	% Yes?
	\draw[line width=1pt, ->, >=stealth] (4.5cm,-11.8cm) -- (4.5cm,-12.7cm) -- (4.5cm,-12.7cm) -- (4.5cm,-13.5cm);
	\draw (4.9cm,-12.0cm) node {\textsf{Yes}};

	% Post-processing.
	\draw (5cm,-13.9cm) node {\textsf{Post-Processing}};
	% Large dashed rectangle.
	\draw[line width=1.2pt, rounded corners, dashed] (0cm,-14.2cm) rectangle (10cm,-17.4cm);
	% Transmission/Reflection coefficient calculations.
	\draw[line width=1pt, rounded corners] (0.5cm,-14.4cm) rectangle (9.5cm,-15.2cm);
	\draw (5cm,-14.8cm) node {\textsf{Compute Transmission and Reflection Coefficients}};
	% Refractive index calculations.
	\draw[line width=1pt, rounded corners] (0.5cm,-15.4cm) rectangle (9.5cm,-16.2cm);
	\draw (5cm,-15.8cm) node {\textsf{Compute Refractive Index of Slab}};
	% Simulation movie.
	\draw[line width=1pt, rounded corners] (0.5cm,-16.4cm) rectangle (9.5cm,-17.2cm);
	\draw (5cm,-16.8cm) node {\textsf{Play Time-Domain Simulation Movie}};
\end{tikzpicture}
\caption{Simulation algorithm}
\label{1D-DNG-Algorithm}
\end{figure}
\section{Simulation Results}
The simulation is run for both lossless and lossy cases with sinusoidal\index{source!sinusoidal}, Gaussian\index{source!Gaussian pulse} and Ricker\index{source!Ricker wavelet} wavelet sources. The slab parameters are set such that at frequency of operation, $f_0$, the permittivity and permeability of slab are both negative and result in a refractive index $n=-1$.
\subsection{Simulation Parameters}
The number of spatial steps is set as $4096$ and simulation is run for $4\times 4096$ time steps. The slab is located between steps $1365$ and $2731$. $\Delta z$ or spatial step is set as $3$ mm and time step, $\Delta t$, is set as $50$ ps. Frequency of operation is $f_0=0.1953125$ GHz and Courant number\index{Courant number ($S_c$)} for this configuration comes out to be $S_c=0.5$. In order to obtain relative permittivity and permeability of $-1$ at required $f_0$, plasma frequencies are set as $\omega^2_{pm}=\omega^2_{pe}=2\times(2\pi f_0)^2$ with $\epsilon_\infty=\mu_\infty=1$. First order absorbing boundary condition (ABC) are applied on fields at end points.
\subsection{Incident and Transmitted Fields}
Simulation with Gaussian pulse\index{source!Gaussian pulse} reveals that low frequency components are reflected at the interface which is confirmed from the transmission and reflection coefficients obtained for the air--slab interface. At $f_0$, the transmission coefficient is $1$ and there are no reflections when a sinusoidal source\index{source!sinusoidal} with $f_0$ is incident on the slab. Under steady--state conditions, transmitted wave inside the slab has negative phase velocity while energy is propagating in $+\hat{z}$ direction as expected.
\subsection{Refractive Index}
Following \cite{DNG-Ehud-Ziol}, the refractive index\index{refractive index ($n$)} was calculated from
\begin{equation}
n_{FDTD} = \dfrac{1}{jk_0(z_1-z_2)}log\left|\dfrac{E_x(\omega,z_2)}{E_x(\omega,z_1)}\right|.
\label{Refractive-Index-FDTD}
\end{equation}
Where, $k_0$ was the wave number\index{wave number ($k_0$)} set as $\omega_0/c$ and the fields were recorded at locations $z_1=1415\Delta z$ and $z_2=1424\Delta z$. For both, Gaussian\index{source!Gaussian pulse} pulse and Ricker\index{source!Ricker wavelet} wavelet, $Re(n)$ was $-1$ at $f_0$ while $Im(n)$ was sufficiently close to $0$.
\begin{figure}[H]
\centering
\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_IncidentFieldGaussian.pdf}
\caption{Gaussian pulse}
\label{1DDNG-IncidentField-Gaussian}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_IncidentFieldRicker.pdf}
\caption{Ricker wavelet}
\label{1DDNG-IncidentField-Ricker}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_IncidentFieldSinusoidal.pdf}
\caption{Sinusoidal wave}
\label{1DDNG-IncidentField-Sinusoidal}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_TransmissionReflectionCoefficient.pdf}
\caption{Transmission and reflection coefficients}
\label{1DDNG-Transmission-Reflection-Coefficient}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_TransmittedField.pdf}
\caption{Transmitted Gaussian pulse}
\label{1DDNG-Transmitted-Gaussian-Pulse}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_TransmittedFieldBeyondSlab.pdf}
\caption{Transmitted Gaussian pulse beyond slab}
\label{1DDNG-Transmitted-Gaussian-Pulse-Beyond-Slab}
\end{figure}
\begin{figure}[H]
\centering
\subfigure{\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_RefractiveIndex.pdf}}
\subfigure{\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_RefractiveIndexZoomed.pdf}}
\caption{Refractive index}
\label{1DDNG-Refractive-Index}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale=0.78, trim=3.5cm 8.7cm 4.5cm 8.75cm, clip]{Figures/FigCh03_1DDNGSteadyStateLossless.pdf}
\caption{Steady-state under lossless conditions}
\label{1DDNG-SteadyState-Lossless}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale=0.78, trim=3.5cm 8.7cm 4.5cm 8.75cm, clip]{Figures/FigCh03_1DDNGSteadyStateLossy.pdf}
\caption{Steady-state under lossy conditions}
\label{1DDNG-SteadyState-Lossy}
\end{figure}
\section{Simulation of 2D DNG Slab}
\subsection{Update Equations}
The most common 2D configurations are $TE^z$\index{polarisation!$TE^z$} or $TM^z$\index{polarisation!$TM^z$} polarisation where the problem space is confined to $xy$--plane. In $TE^z$, electric field components are transverse to $z$--axis and vice versa. For 2D simulation of DNG slab, $TM^z$ polarisation is assumed. The field components of interest would be $E_z$, $H_x$ and $H_y$. From equations \ref{Faraday's Law} and \ref{Ampere's Law},
\begin{equation}
\nabla \times \textbf{E} = -\dfrac{\partial \textbf{B}}{\partial t} = \left| \begin{array}{ccc} \hat{x} & \hat{y} & \hat{z} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ 0 & 0 & E_z \end{array} \right| = \dfrac{\partial E_z}{\partial y} \hat{x} - \dfrac{\partial E_z}{\partial x} \hat{y}
\label{eq:Faraday'sLawExpansionTMz}
\end{equation}
and
\begin{equation}
\nabla \times \textbf{H} = \dfrac{\partial \textbf{D}}{\partial t} = \left| \begin{array}{ccc} \hat{x} & \hat{y} & \hat{z} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ H_x & H_y & 0 \end{array} \right| = \dfrac{\partial H_y}{\partial x} \hat{z} - \dfrac{\partial H_x}{\partial y} \hat{z}.
\label{eq:Ampere'sLawExpansionTMz}
\end{equation}
The field components are given by scalar equations
\begin{equation}
\dfrac{\partial D_z}{\partial t} = \dfrac{\partial H_y}{\partial x} - \dfrac{\partial H_x}{\partial y},
\label{eq:partialofDzbyt}
\end{equation}
\begin{equation}
\dfrac{\partial B_x}{\partial t} = -\dfrac{\partial E_z}{\partial y}
\label{eq:partialofBxbyt}
\end{equation}
and
\begin{equation}
\dfrac{\partial B_y}{\partial t} = \dfrac{\partial E_z}{\partial x}.
\label{eq:partialofBybyt}
\end{equation}
FDTD update equations are obtained by discretising equations\index{FDTD!update equations} \ref{eq:partialofDzbyt}, \ref{eq:partialofBxbyt} and \ref{eq:partialofBybyt}, given by
\begin{equation}
\begin{split}
D^{n+1}_z \left[i,j\right]=D^{n}_z \left[i,j\right]+\dfrac{\Delta t}{\Delta}&\left(H^{n+\frac{1}{2}}_y\left[i+\frac{1}{2},j\right]-H^{n+\frac{1}{2}}_y \left[i-\frac{1}{2},j\right]\right.\\
&\left.-H^{n+\frac{1}{2}}_x \left[i,j+\frac{1}{2}\right]+H^{n+\frac{1}{2}}_x \left[i,j-\frac{1}{2}\right]\right),
\end{split}
\label{eq:Dz-2D-FDTD-TMz}
\end{equation}
\begin{equation}
B^{n+\frac{1}{2}}_x \left[i,j+\frac{1}{2}\right]=B^{n-\frac{1}{2}}_x \left[i,j+\frac{1}{2}\right] + \dfrac{\Delta t}{\Delta} \left(-E^{n}_z \left[i,j+1\right] + E^{n}_z \left[i,j\right] \right)
\label{eq:Bx-2D-FDTD-TMz}
\end{equation}
and
\begin{equation}
B^{n+\frac{1}{2}}_y \left[i+\frac{1}{2},j\right]=B^{n-\frac{1}{2}}_y \left[i+\frac{1}{2},j\right] + \dfrac{\Delta t}{\Delta} \left( E^{n}_z \left[i+1,j\right] - E^{n}_z \left[i,j\right] \right).
\label{eq:By-2D-FDTD-TMz}
\end{equation}
It is important to note here that DNG slab medium is assumed anisotropic. The auxiliary update equations\index{FDTD!auxiliary update equation} (\ref{2nd-order-B-H-final-form} and \ref{2nd-order-D-E-final-form}) to obtain electric and magnetic fields, $E_v$ and $H_v$ ($v\in x,y,z$), from corresponding flux densities, $D_v$ and $B_v$, remain unchanged.
\subsection{Simulation Geometry}
The DNG slab interface is perpendicular to incident plane wave\index{plane wave} propagating in $+y$ direction. Periodic boundary conditions (PBC)\index{boundary condition!periodic (PBC)} are applied at $x=0$ and $x=x_{max}$. In $y$ direction the grid is terminated at both ends by perfectly matched layer (PML)\index{boundary condition!PML} to absorb any incoming waves. The solution geometry is depicted in figure \ref{fig:2D-DNG-Geometry}. The arrangement of field nodes is shown in figure \ref{fig:2D-DNG-Field-Nodes}.  Magnetic field components are updated first and then used to update electric field.

The solution space is bounded in $y$ direction at both ends by $H_x$ which acts as the boundary for PML. $H_x$ at these end points is set to 0 so that any out--going fields are reflected back into PML. Essentially, the $j$ size of $H_x$ arrays is one greater than $H_y$ and $E_z$ arrays. This arrangement results in following update equations that can be directly programmed, taking into account discrete indices for field storage.
\begin{equation}
\begin{split}
D^{n+1}_z \left[i,j\right]=D^{n}_z \left[i,j\right]+\dfrac{\Delta t}{\Delta}&\left(H^{n+1}_y\left[i,j\right]-H^{n+1}_y \left[i-1,j\right]\right.\\
&\left.-H^{n+1}_x \left[i,j+1\right]+H^{n+1}_x \left[i,j\right]\right).
\end{split}
\label{eq:Dz-2D-FDTD-TMz-Corrected}
\end{equation}
\begin{equation}
B^{n+1}_x \left[i,j\right]=B^n_x \left[i,j\right] + \dfrac{\Delta t}{\Delta} \left(-E^{n}_z \left[i,j\right] + E^{n}_z \left[i,j-1\right] \right).
\label{eq:Bx-2D-FDTD-TMz-Corrected}
\end{equation}
\begin{equation}
B^{n+1}_y \left[i,j\right]=B^n_y \left[i,j\right] + \dfrac{\Delta t}{\Delta} \left( E^{n}_z \left[i+1,j\right] - E^{n}_z \left[i,j\right] \right).
\label{eq:By-2D-FDTD-TMz-Corrected}
\end{equation}
\subsection{Simulation Parameters and Results}
The solution space is 512 cells in both $x$ and $y$ directions without taking into account the width of PML, which is 50 cells wide. The plane wave source is located 10 cells from the lower PML layer. The spatial and temporal steps are set as $\Delta x=\Delta y=\Delta=3$ mm and $\Delta t=50$ ps, respectively; with a Courant number of 0.5\index{Courant number ($S_c$)}. Frequency of operation is $f_0=1.5625$ GHz. The DNG slab parameters are same as in the case of 1D DNG simulation.

Figure \ref{fig:2DDNG-Refractive-Index} shows the refractive index\index{refractive index ($n$)} obtained for a range of frequencies with a Gaussian pulse excitation. Again, real part of refractive index at $f_0$ is close to -1, whereas, imaginary part is close to 0. Figure \ref{fig:2DDNG-SteadyState-Lossy} shows sinusoidal steady--state\index{steady--state} under lossy conditions where the source is located towards left and incident plane wave is moving rightwards. Under lossless conditions the slab can act as a perfect lens\index{perfect lens} (figure \ref{fig:2DDNG-Slab-lens}).
\begin{figure}[H]
\centering
\begin{tikzpicture}[xscale=0.9,yscale=0.9]
	\newcommand{\LeftX}{0cm}
	\newcommand{\MidX}{4cm}
	\newcommand{\RightX}{8cm}
	\newcommand{\PMLw}{0.8cm}
	\newcommand{\DomainY}{6.0cm}
	\newcommand{\SlabStartY}{2.0cm+\PMLw}
	\newcommand{\SlabEndY}{4.0cm+\PMLw}
	% x-axis.
	\draw[->, >=stealth] (\LeftX,0cm) -- (\RightX+0.5cm,0cm);
	\coordinate [label=right:$x$] (x-axis) at (\RightX+0.6cm,0cm);
	% y-axis.
	\draw[->, >=stealth] (\LeftX,0cm) -- (\LeftX,2*\PMLw+\DomainY+0.5cm);
	\coordinate [label=above:$y$] (y-axis) at (\LeftX,2*\PMLw+\DomainY+0.6cm);
	% PBCs.
	\coordinate [label=right:PBC] (PBCright) at (\RightX,\PMLw+3.0cm);
	\coordinate [label=left:PBC] (PBCleft) at (\LeftX,\PMLw+3.0cm);
	% Lower PML Region.
	\draw[fill=gray!40!white] (\LeftX,0cm) rectangle (\RightX,\PMLw);
	\coordinate [label=center:\textsf{Lower PML Region}] (LowerPML) at (\MidX,0.4cm);
	% Solution Region.
	\draw (\LeftX,\PMLw) rectangle (\RightX,\PMLw+\DomainY);
	% Upper PML Region.
	\draw[fill=gray!40!white] (\LeftX,\PMLw+\DomainY) rectangle (\RightX,2*\PMLw+\DomainY);
	\coordinate [label=center:\textsf{Upper PML Region}] (UpperPML) at (\MidX,\PMLw+\DomainY+0.4cm);
	% Slab.
	\draw[fill=gray!10!white] (\LeftX,\SlabStartY) rectangle (\RightX,\SlabEndY);
	\coordinate [label=center:\textsf{DNG Slab}] (DNGSlab) at (\MidX,\PMLw+3.0cm);
	% Plane wave.
	\draw (\LeftX+2.0cm,\PMLw+0.2cm) -- (\RightX-2.0cm, \PMLw+0.2cm);
	\draw (\LeftX+2.0cm,\PMLw+0.3cm) -- (\RightX-2.0cm, \PMLw+0.3cm);
	\draw (\LeftX+2.0cm,\PMLw+0.4cm) -- (\RightX-2.0cm, \PMLw+0.4cm);
	\draw[line width=1.1pt, ->, >=stealth] (\MidX,\PMLw+0.6cm) -- (\MidX,\PMLw+1.6cm);
	\coordinate [label=right:\textsf{Plane Wave}] (PlaneWave) at (\MidX+0.5cm,\PMLw+0.8cm);
\end{tikzpicture}
\caption{Simulation geometry}
\label{fig:2D-DNG-Geometry}
\end{figure}
\begin{figure}[H]
\centering
\vspace{-0.5cm}
\begin{tikzpicture}
	\newcommand{\LeftX}{0cm}
	\newcommand{\RightX}{5cm}
	\newcommand{\LowY}{0cm}
	\newcommand{\HighY}{5.5cm}
	% x-axis.
	\draw[->, >=stealth] (\LeftX,\LowY) -- (\RightX,\LowY);
	\coordinate [label=right:$x$] (x-axis) at (\RightX+0.1cm,\LowY);
	% y-axis.
	\draw[->, >=stealth] (\LeftX,\LowY) -- (\LeftX,\HighY);
	\coordinate [label=above:$y$] (y-axis) at (\LeftX,\HighY+0.1cm);
	% Drawing vertical grid lines.
	\foreach \x in {1cm,2cm,3cm,4cm}
		\draw (\x,\LowY-0.1cm) -- (\x,\HighY-0.25cm); % Solid lines at +1 intervals.
	\foreach \x in {0.5cm,1.5cm,2.5cm,3.5cm,4.5cm}
		\draw[dashed] (\x,\LowY) -- (\x,\HighY-0.25cm); % Dashed lines at 1/2 intervals.
	% Drawing horizontal grid lines.
	\foreach \y in {1cm,2cm,3cm,4cm,5cm}
		\draw (\LeftX-0.1cm,\y) -- (\RightX-0.25cm,\y); % Solid lines at +1 intervals.
	\foreach \y in {0.5cm,1.5cm,2.5cm,3.5cm,4.5cm}
		\draw[dashed] (\LeftX,\y) -- (\RightX-0.25cm,\y); % Dashed lines at 1/2 intervals.
	% Drawing nodes.
	\foreach \x in {0cm,1cm,2cm,3cm,4cm}
		\foreach \y in {0.5cm,1.5cm,2.5cm,3.5cm,4.5cm}
			\filldraw (\x,\y) circle (0.08cm); % Ez nodes.
	\foreach \x in {0.5cm,1.5cm,2.5cm,3.5cm,4.5cm}
		\foreach \y in {0.5cm,1.5cm,2.5cm,3.5cm,4.5cm}
			\node[fill=black,regular polygon, regular polygon sides=4,inner sep=0.05cm] at (\x,\y) {}; % Hy nodes.
	\foreach \x in {0cm,1cm,2cm,3cm,4cm}
		\foreach \y in {0cm,1cm,2cm,3cm,4cm,5cm}
			\node[fill=black,regular polygon, regular polygon sides=3,inner sep=0.04cm] at (\x,\y) {}; % Hx nodes.
	% +1 Text for x-axis.
	\foreach \x/\t in {0cm/0,1cm/1,2cm/2,3cm/3,4cm/4}
		\coordinate [label=below:$\t$] (\t) at (\x,\LowY-0.2cm);
	% 1/2 Text for x-axis.
	\foreach \x/\t in {0.5cm/ ,1.5cm/1,2.5cm/2,3.5cm/3,4.5cm/4}
		\coordinate [label=below:{\tiny $\t\frac{1}{2}$}] (\t) at (\x,\LowY-0.2cm);
	% +1 Text for y-axis.
	\foreach \y/\t in {0cm/0,1cm/1,2cm/2,3cm/3,4cm/4,5cm/5}
		\coordinate [label=left:$\t$] (\t) at (\LeftX-0.2cm,\y);
	% 1/2 Text for y-axis.
	\foreach \y/\t in {0.5cm/ ,1.5cm/1,2.5cm/2,3.5cm/3,4.5cm/4}
		\coordinate [label=left:{\tiny $\t\frac{1}{2}$}] (\t) at (\LeftX-0.2cm,\y);
	% Legend
	\filldraw (\LeftX+0cm,\LowY-1cm) circle (0.08cm); % Ez.
	\coordinate [label=center:$E_z$] (EzLegend) at (\LeftX+0.4cm,\LowY-1cm);
	\node[fill=black,regular polygon, regular polygon sides=3,inner sep=0.04cm] at (\LeftX+1cm,\LowY-1cm) {}; % Hx.
	\coordinate [label=center:$H_x$] (HxLegend) at (\LeftX+1.4cm,\LowY-1cm);
	\node[fill=black,regular polygon, regular polygon sides=4,inner sep=0.05cm] at (\LeftX+2cm,\LowY-1cm) {}; % Hy.
	\coordinate [label=center:$H_y$] (HxLegend) at (\LeftX+2.4cm,\LowY-1cm);
\end{tikzpicture}
\caption{Arrangement of field nodes}
\label{fig:2D-DNG-Field-Nodes}
\end{figure}
\begin{figure}[H]
\centering
\subfigure{\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_2DDNGRefractiveIndex.pdf}}
\subfigure{\includegraphics[scale=0.8, trim=3.5cm 8.7cm 4.5cm 8.85cm, clip]{Figures/FigCh03_2DDNGRefractiveIndexZoomed.pdf}}
\caption{Refractive index of 2D DNG slab}
\label{fig:2DDNG-Refractive-Index}
\end{figure}
\begin{figure}[H]
\centering
\subfigure{\includegraphics[scale=0.5, trim=1.0cm 0.05cm 1.75cm 0cm, clip]{Figures/FigCh03_2DDNGSteadyStateLossy.png}}
\subfigure{\includegraphics[scale=0.5, trim=1.0cm 0.05cm 1.75cm 0cm, clip]{Figures/FigCh03_2DDNGSteadyStateLossyOverhead.png}}
\caption{Steady-state under lossy conditions for 2D DNG slab}
\label{fig:2DDNG-SteadyState-Lossy}
\end{figure}
\begin{figure}[H]
\centering
\subfigure{\includegraphics[scale=0.5, trim=1.0cm 0.05cm 1.75cm 0cm, clip]{Figures/FigCh03_2DDNGSlabLensSideView.png}}
\subfigure{\includegraphics[scale=0.5, trim=1.0cm 0.05cm 1.75cm 0cm, clip]{Figures/FigCh03_2DDNGSlabLens.png}}
\caption{DNG slab as lens: source is located towards left}
\label{fig:2DDNG-Slab-lens}
\end{figure}
